Zadanie 1

local f16 = Float16(1/3) ::Float16

Float16(0.3333)

local f32 = Float32(1/3) ::Float32

0.33333334f0

local f64 = Float64(1/3) ::Float64

0.3333333333333333

Porównanie reprezentacji bitowej liczby 1/3 dla Float16, Float32, Float64

bitstring(f16)

"0011010101010101"

bitstring(f32)

"00111110101010101010101010101011"

bitstring(f64)

"0011111111010101010101010101010101010101010101010101010101010101"

Float32(1/3) == Float16(1/3)

false

Float32(1/3) == Float64(1/3)

false

Float16(1/3) == Float64(1/3)

false

Ta sama liczba zmiennoprzecinkowa w różnych reprezentacjach zmiennoprzecinkowych nie jest sobie równa.

Liczba zainicjowana jako Float16 a potem rzutowana na Float64

convertedTof64 = convert(Float64, f16)

0.333251953125

f64 == convertedTof64

false

Uwaga Float16 jest równy Floatowi16 rzutowanemu na Float64

f16 == convertedTof64

true

f32 == convertedTof64

false

Zadanie 2

Maszynowe epsilon mówi, jaki jest odstęp między kolejnymi liczbami zmiennoprzecinkowymi reprezentowanymi w komputerze

arguments = Float32[]
values = Float32[]
_start = 1.0
_end = 1000000.0
range = _start : _end
for i=range
    push!(arguments, i)
    push!(values, eps(i))
end

using Plots
plot(arguments, values)

Zauważamy że im dalej od zera tym większe są odległości między liczbami zmiennoprzecinkowymi.

Zadanie 3

#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>

int main(){
    float a = 1e-33;
    while(a > 0){
        a /= 2.0;
        printf (" f="); gsl_ieee_printf_float(&a);
        printf("\n");
    }
    return 0;
}

 Fragment wyniku:

 f= 1.01001100010011101001100*2^-122
 f= 1.01001100010011101001100*2^-123
 f= 1.01001100010011101001100*2^-124
 f= 1.01001100010011101001100*2^-125
 f= 1.01001100010011101001100*2^-126
 f= 0.10100110001001110100110*2^-126  // Tutaj mantysa przestaje być znormalizowana
 f= 0.01010011000100111010011*2^-126
 f= 0.00101001100010011101010*2^-126
 f= 0.00010100110001001110101*2^-126
 f= 0.00001010011000100111010*2^-126
 f= 0.00000101001100010011101*2^-126
 f= 0.00000010100110001001110*2^-126

Domyślny bit w mantysie przestaje mieć wartość 1, z tego wnioskujemy że liczba przestaje być znormalizowana

Zadanie 4

Gdy liczymy miejsce zerowe równania kwadratowego ax^2 + bx + c dla a bliskiego 0, ze wzoru x = (-b + sqrt(b^2 - 4ac)) / 2a <1> , dostajemy stosunkowo duże błedy wyniku. Możemy przekształcić wzór <1> za pomocą równania a-b = (a^2 - b^2)/ (a + b), więc x = (b^2 - 4ac - b^2)/(2a (b + sqrt(b^2 - 4 a c)) => x = -2c / (b + sqrt(b^2 - 4 a c)) wtedy równanie już nie jest tak czułe na zaburzenie a.

unstableQuadratic(a, b, c) =  (-b + sqrt(b^2 - 4*a*c)) / 2a

unstableQuadratic (generic function with 1 method)

stableQuadratic(a, b, c) = -2*c / (b + sqrt(b^2 - 4 * a * c))

stableQuadratic (generic function with 1 method)

valuesUnstable = Float32[]
valuesStable = Float32[]
arguments = Float32[]
a = 5.0 * 10 ^ -15
b = 1
c = -2
while a > 0
    push!(valuesUnstable, unstableQuadratic(a, b, c))
    push!(valuesStable, stableQuadratic(a, b, c))
    push!(arguments, a)
    a = a / 2
end

Dla podanych paramterów (a = 5.0 * 10 ^ -15, b = 1, c = -2) ax^2 + bx + c => ax^2 + x - 2, dla a bliskiego 0, x powinno być równe 2. Dla algorytmu niestabilnego tak nie jest.

Wersja niestabilna

using Plots
plot(arguments, valuesUnstable, label="Unstable")

Wersja stabilna w porównaniu do niestabilnej

plot!(arguments, valuesStable, label="Stable")

Dla a bliskiego 0, algorytm niestabilny daje wynik 0, który jest daleki od poprawnego 2. Powodem tego błędu jest dzielenie przez bardzo małą liczbe.